********************************************************************************
*
*	Understanding the Impact of US Military Service on Support for Insurrection
*   Journal of Conflict Resolution
*   Replication
*
*   Robert A. Pape
*	Keven G. Ruby   
*	Kyle D. Larson
*   Kentaro Nakamura
*
*	Replicate supplement figures and tables (except for table a7)
*   
*   v2024-05-29
*
********************************************************************************


********************************************************************************
* 
* Setup
*
********************************************************************************

// Stata version = StataNow 18.5
version 18.5

// These user-installed packages are needed for this do-file to run properly 
ssc install catplot, replace 
net install sensemakr,  replace force from("https://raw.githubusercontent.com/resonance1/sensemakr-stata/master/")
ssc install catplot, replace 
ssc install coefplot, replace 
ssc install schemepack, replace
ssc install blindschemes, replace
ssc install heatplot, replace


// Scheme used for figures
set scheme white_tableau

* Load dataset
use vetstudy_jcr_use.dta, clear


// Set Survey Weighting for comparisons between unweighted and weighted regressions
svyset vpsu32, strata(vstrat32) weight(weight_all) vce(linearized) singleunit(center) 

/* Establish List of Pre-Treatment Controls for All Subsequent Regressions 
local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white" // only binary controls
local controls_c "b_male age b_christian religiosity childlaw b_republican b_independent b_white" // includes continous controls
local mediators "1.b_conspiracy 1.b_polarization"
local servicevars "1.b_combattrauma 1.b_discharged 1.b_servedgwot 1.b_singletour 1.b_combatdeployment"
local postservicevars "1.b_unappreciated 1.b_notproud 1.b_reintegration 1.b_disillusioned 1.b_ptsd"		
*/

********************************************************************************
*
* Table A.1. 
* Descriptive Statistics for Dependent Variable and Covariates
*
********************************************************************************

dtable insurrection_index `controls' 1.b_authoritarian `mediators' 1.b_econconcern 1.b_hardship `servicevars' `postservicevars', ///
by(veteran, nototals) export(tab_a1.docx, replace)


********************************************************************************
*
* 	Figure A.1.
*   Correlation matrix for Political Conspiracy Index and components (unweighted)
*
********************************************************************************

pwcorr conspiracy_qanon conspiracy_jailpatriots conspiracy_demreplacement conspiracy if veteran==1
matrix corrmatrix = r(C)
heatplot corrmatrix, values(format(%9.2f) size(vsmall)) color(hcl diverging, intensity(.6)) xlabel(,angle(90)) legend(off) aspect(1) title(Veterans) lower cuts(-1.05(.1)1.05) name(vets, replace) xsize(2) ysize(2)

pwcorr conspiracy_qanon conspiracy_jailpatriots conspiracy_demreplacement conspiracy if veteran==0
matrix corrmatrix = r(C)
heatplot corrmatrix, values(format(%9.2f) size(vsmall)) color(hcl diverging, intensity(.6)) xlabel(,angle(90)) legend(off) aspect(1) title(Non-Veterans) lower cuts(-1.05(.1)1.05) name(nvets, replace) xsize(2) ysize(2)

graph combine vets nvets
graph export fig_a1.png, replace width(3000) 


********************************************************************************
*
* 	Figure A.2.
*   Unweighted distribution of Political Conspiracy Index (veteran vs non-veteran samples)
*
********************************************************************************

gen rconspiracy = round(conspiracy,.1)

catplot veteran, over(rconspiracy, label(labsize(medium))) 					///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 					///
	ytitle("Percent of Sample", size(medium)) 								///
	b1title("Far-Right Conspiracy Index") 			        				///
	bargap(3) scheme(white_tableau) percent(veteran) 						///
	vertical blabel(bar, format(%4.0f) size(medium))                     	///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a2.png, replace width(3000)


********************************************************************************
*
* 	Figure A.3.
*   Unweighted distribution of Believing Far-Right Conspiracy Index (veteran vs non-veteran samples)
*
********************************************************************************

catplot veteran, over(b_conspiracy) 						///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 							///
	l1title("Believe Far-Right Conspiracy" "(Far-Righ Conspiracy Index >= 3.5)") 			///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	horizontal blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a3.png, replace width(3000)

********************************************************************************
*
* 	Figure A.4.
*   Unweighted distribution of Relative Support for Trump (veteran vs non-veteran samples)
*
********************************************************************************
 
catplot veteran, over(polarization, label(labsize(medium))) 						    ///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 							///
	b1title("Political Polarization Index") 			        ///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	vertical blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a4.png, replace width(3000)


********************************************************************************
*
* Figure A.5. Unweighted distribution of Relative Support for Trump (veteran vs non-veteran samples)
*
********************************************************************************

catplot veteran, over(b_polarization) 						///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 							///
	l1title("Pro-Trump Polarization (Polarization Index >= 3)") 			///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	horizontal blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a5.png, replace width(3000)	


********************************************************************************
*
* Figure A.6. Correlation matrix for Combat Experience and components for veteran and non-veteran samples (unweighted)
*
********************************************************************************

pwcorr combat mortalfear sawfriendkilled sawenemykilled sawcivskilled if veteran==1
matrix corrmatrix = r(C)
heatplot corrmatrix, values(format(%9.2f) size(vsmall)) color(hcl diverging, intensity(.6)) xlabel(,angle(90)) legend(off) aspect(1) title(Veterans) lower cuts(-1.05(.1)1.05) name(vets, replace) xsize(2) ysize(2)
graph export fig_a6.png, replace width(3000)


********************************************************************************
*
* Figure A.7. Unweighted distribution of Combat Experience (veteran sample only)
*
********************************************************************************

catplot combattrauma if veteran==1, percent 						///
	bar(1, color(green*.7))  	///
	ytitle("Percent of Sample", size(medium)) 							///
	b1title("Combat Trauma Index") 			        ///
	bargap(3) scheme(white_tableau)  		///
	vertical blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a7.png, replace width(3000)


********************************************************************************
*
* Figure A.8. Unweighted distribution of Experienced Combat Trauma (Combat Trauma Index > 1) (veteran sample only)
*
********************************************************************************

catplot b_combattrauma if veteran==1, percent 						///
	bar(1, color(green*.7))	///	///
	ytitle("Percent of Sample", size(medium)) 							///
	l1title("Experienced Combat Trauma (Combat Index > 1)") 			///
	bargap(3) scheme(white_tableau) 		///
	horizontal blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a8.png, replace width(3000)

********************************************************************************
*
* Figure A.9. Unweighted Distribution of Served During GWoT (veteran sample only)
*
********************************************************************************
catplot b_servedgwot if veteran==1, percent 						///
	bar(1, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 							///
	l1title("Served in GWOT (post 2001)") 			///
	bargap(3) scheme(white_tableau) 		///
	horizontal blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a9.png, replace width(3000)

********************************************************************************
*
* Figure A.10. Correlation matrix for Religiosity and components for veteran and non-veteran samples (unweighted)
*
********************************************************************************

pwcorr relig_pray relig_services relig_life religiosity if veteran==1
matrix corrmatrix = r(C)
heatplot corrmatrix, values(format(%9.2f) size(vsmall)) color(hcl diverging, intensity(.6)) xlabel(,angle(90)) legend(off) aspect(1) title(Veterans) lower cuts(-1.05(.1)1.05) name(vets, replace) xsize(2) ysize(2)

pwcorr relig_pray relig_services relig_life religiosity if veteran==0
matrix corrmatrix = r(C)
heatplot corrmatrix, values(format(%9.2f) size(vsmall)) color(hcl diverging, intensity(.6)) xlabel(,angle(90)) legend(off) aspect(1) title(Non-Veterans) lower cuts(-1.05(.1)1.05) name(nvets, replace) xsize(2) ysize(2)
graph combine vets nvets
graph export fig_a10.png, replace width(3000)

********************************************************************************
*
* Figure A.11. Unweighted distribution of Religiosity (veteran vs non-veteran samples)
*
********************************************************************************
replace religiosity = round(religiosity,0.1)

catplot veteran, over(religiosity) 						    ///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 							///
	b1title("Religiosity Index") 			        ///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	vertical blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a10.png, replace width(3000)

********************************************************************************
*
* Figure A.12. Unweighted distribution of Religious (Religiosity Index > 4) (veteran vs non-veteran samples)
*
********************************************************************************

catplot veteran, over(b_religious) 						///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 							///
	l1title("Religious (Religiosity Index > 4)") 			///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	horizontal blabel(bar, format(%4.0f) size(medium))                     ///
	xsize(6) ysize(2) legend(size(medium)) ylab(,labsize(medium))
graph export fig_a12.png, replace width(3000)


estimates clear 

********************************************************************************
*
* Figure A.13. Unweighted distribution the Insurrectionist Sentiments Index and three components (veteran vs non-veteran samples)
*
********************************************************************************
 

catplot veteran, over(insurrection_index) 				    ///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 				///
	b1title("Insurrection Sentiments Index") 			    ///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	vertical blabel(bar, format(%4.0f) size(medium))        ///
	xsize(6) ysize(2) legend(size(medium) pos(3) row(2)) ylab(,labsize(medium)) ///
    name(gindex, replace)

catplot veteran, over(force, nolab) 				    	///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 	///
	ytitle("Percent of Sample", size(medium)) 				///
	b1title("Force Justified to Restore Trump") 			///
	bargap(3) scheme(white_tableau) percent(veteran) 		///
	vertical blabel(bar, format(%4.0f) size(medium))       	///
	xsize(6) ysize(2) legend(off) ylab(,labsize(medium)) 	///
    name(gforce, replace)	
	
catplot veteran, over(force_personal, nolab) 					///
	asyvar bar(1, color(blue*.7)) bar(2, color(green*.7)) 		///
	ytitle("Percent of Sample", size(medium)) 					///
	b1title("I Would Personally Use Force to Restore Trump") 	///
	bargap(3) scheme(white_tableau) percent(veteran) 			///
	vertical blabel(bar, format(%4.0f) size(medium)) 			///
	xsize(6) ysize(2) legend(off) ylab(,labsize(medium)) 		///
    name(gpersonal, replace)	

graph combine gforce gpersonal, ycommon name(gtoprow, replace)
graph combine gtoprow gindex, row(2) ycommon 

graph export fig_a13.png, replace width(3000)


********************************************************************************
*
* Table A.2. Effect of Veteran on Index Components
*
********************************************************************************
local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white"

regress insurrection_index veteran `controls', robust
estimates store tab6
 
 regress insurrection_index veteran `controls' b_authoritarian, robust
estimates store tab6_b

 regress insurrection_index veteran `controls' b_nocollege , robust
estimates store tab6_e

 
// Across All Specifications of DV
foreach dv of varlist force force_personal {

	regress `dv' veteran `controls', robust
	estimates store tab6_`dv'

	regress `dv' veteran `controls' b_authoritarian, robust
	estimates store tab6_`dv'_b	

	regress `dv' veteran `controls' b_nocollege, robust
	estimates store tab6_`dv'_e	
	
}

esttab tab6 tab6_b tab6_e tab6_force tab6_force_b tab6_force_e tab6_force_personal tab6_force_personal_b tab6_force_personal_e using tab_a2.rtf ///
, replace se compress r2 label ///
mtitles("Index" "w/ Authoritarian" "w/ Education" "Force" "w/ Authoritarian" "w/ Education" "Personal" "w/ Authoritarian" "w/ Education") ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f)  ///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: Unweighted OLS with robust standard errors. Standard errors in parentheses.")

********************************************************************************
*
* Table A.3. Effect of Veteran using Engaged and Weighting
*
********************************************************************************

/* Inattentive Variables

	Inconsistent Answers 	Answered two nearly-identical questions oppositely
	Duration 				More than one standard deviation faster than mean
	

*/

gen inconsistent = 0
	/// doubt 2020 election results & strongly disagree 2020 election stolen
    replace inconsistent = 1 if doubt == 1 & stolen == 1 	 
	/// doubt 2020 election results & somewhat disagree 2020 election stolen
    replace inconsistent = 1 if doubt == 1 & stolen == 2 
	/// trust 2020 election results & somewhat agree 2020 election stolen
    replace inconsistent = 1 if doubt == 2 & stolen == 4 	 
	/// trust 2020 election results & strongly agree 2020 election stolen
    replace inconsistent = 1 if doubt == 2 & stolen == 5 	 
								
egen zduration = std(duration), by(veteran)		// Standard Deviations on Duration separately for vets and non-vets
gen rzduration = round(zduration)				// Rounded Standardized Duration
	replace rzduration = -1.5 if zduration <-1	// More than one SD faster relative to group mean (vets/non-vets)
	replace rzduration = 3 if rzduration >3 	// More than three SD slower 

gen speeder = 0 							
	replace speeder = 1 if rzduration == -1.5	// Speeder defined as > 1 SD faster
	replace speeder = . if missing(duration)

gen disengaged = 0 
	replace disengaged = 1 if speeder == 1 | inconsistent == 1


local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white"

// Main - Unweighted
regress insurrection_index veteran `controls', robust
estimates store tab6

// Engaged Only - Unweighted
regress insurrection_index veteran `controls' ///
	if disengaged == 0 , robust
estimates store tab6_b

// Main - Weighted
svy: regress insurrection_index veteran `controls'
estimates store tab6_w

// Engaged - Weighted
svy: regress insurrection_index veteran `controls' ///
	if disengaged == 0 
estimates store tab6_wb


esttab tab6 tab6_b tab6_w tab6_wb using tab_a3.rtf ///
, replace se compress r2 label ///
mtitles("Main" "Engaged" "Main" "Engaged") ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
coeflabels(_cons "Constant" veteran "Veteran") ///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: OLS, robust standard errors in parentheses")
  
******************************************************************************** 
*
* Table A.4. Effect of Veteran using Less Restrictive Specifications of IVs 
*
********************************************************************************
local controls_c "b_male age b_christian religiosity childlaw b_republican b_independent b_white" // includes continous controls

// Main - Unweighted
regress insurrection_index veteran `controls_c', robust
estimates store tab6c

// Engaged Only - Unweighted
regress insurrection_index veteran `controls_c' ///
	if disengaged == 0 , robust
estimates store tab6c_b

// Main - Weighted
svy: regress insurrection_index veteran `controls_c'
estimates store tab6c_w

// Engaged - Weighted
svy: regress insurrection_index veteran `controls_c' ///
	if disengaged == 0 
estimates store tab6c_wb


esttab tab6c tab6c_b tab6c_w tab6c_wb using tab_a4.rtf ///
, replace se compress r2 label ///
mtitles("Main" "Engaged" "Main" "Engaged") ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
coeflabels(_cons "Constant" veteran "Veteran") ///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: OLS, robust standard errors in parentheses")

********************************************************************************  
*
* Table A.5. Effect of Veterancy on Having High Insurrectionist Sentiments  
*
********************************************************************************
local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white"

// Main - Unweighted
logit b_insurrect veteran `controls', robust or
estimates store tab6b

// Engaged Only - Unweighted
logit b_insurrect veteran `controls' ///
	if disengaged == 0 , robust or
estimates store tab6b_b

// Main - Weighted
svy: logit b_insurrect veteran `controls', or
estimates store tab6b_w

// Engaged - Weighted
svy: logit b_insurrect veteran `controls' ///
	if disengaged == 0, or 
estimates store tab6b_wb


esttab tab6b tab6b_b tab6b_w tab6b_wb using tab_a5.rtf ///
, replace compress r2 label ///
mtitles("Main" "Engaged" "Main" "Engaged") ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
coeflabels(_cons "Constant" veteran "Veteran") ///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: logistic regression, coefficients reported as odds ratios. Robust standard errors in parentheses")

********************************************************************************  
*
* Table A.6. Effect of Veterancy on Having High Insurrectionist Sentiments (Less Restrictive IV Specifications)
*
********************************************************************************

local controls_c "b_male age b_christian religiosity childlaw b_republican b_independent b_white" // includes continous controls

// Main - Unweighted
logit b_insurrect veteran `controls_c', robust or
estimates store tab6bc

// Engaged Only - Unweighted
logit b_insurrect veteran `controls_c'  ///
	if disengaged == 0 , robust or
estimates store tab6bc_b

// Main - Weighted
svy: logit b_insurrect veteran `controls_c', or
estimates store tab6bc_w

// Engaged - Weighted
svy: logit b_insurrect veteran `controls_c' ///
	if disengaged == 0, or 
estimates store tab6bc_wb


esttab tab6bc tab6bc_b tab6bc_w tab6bc_wb using tab_a6.rtf ///
, replace se compress label ///
mtitles("Main" "Engaged" "Main" "Engaged") ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
coeflabels(_cons "Constant" veteran "Veteran") ///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: logistic regression, coefficients reported as odds ratios. Robust standard errors in parentheses")
  

  
********************************************************************************  
*
* Figure A.14. Sensitivity Analysis with Less Restrictive Specifications of the IVs. 
*
********************************************************************************
local controls_c "b_male age b_christian religiosity childlaw b_republican b_independent b_white" // includes continous controls

regress insurrection_index veteran `controls_c', robust 

 
sensemakr 	insurrection_index 									 /// DV
			veteran 											 /// Primary IV
			`controls_c' 									     ///
	, treat(veteran)  											 /// Treatment
	benchmark													 /// Benchmarks
		(childlaw age b_christian b_republican) ///
	contourplot kd(4) clim(0 0.3)

graph display, xsize(1.5) ysize(1.5) 		  
graph export fig_a14.png, replace width(3000)


drop if missing(caseid) // This is a consequence of the aftermath of the 
						// Sensemaker package, which creates additional records
						// that are now dropped so we can resume analysis 

********************************************************************************
*
* Figure A.15. Average Controlled Direct Effect of Veteran on Insurrectionist Sentiments by Mediator After Dropping Disengaged
*
********************************************************************************

// Mediation Analysis, Dropping Inattentive Respondents 
preserve  
	
	drop if disengaged==1
	capture gen ytilde=.
	
/* Citation
Method using:

Acharya, Avidit, Matthew Blackwell, and Maya Sen. 2016. "Explaining Causal 
Findings Without Bias: Detecting and Assessing Direct Effects." American 
Political Science Review 110 (3): 512–29. https://doi.org/10/f9cr9k.
*/

local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white"

* Define Variables Used in the Mediation Analysis, by Category *****************
local y  "insurrection_index" 			// Dependent Variable 
local a  "veteran"						// Treatment 
local z  "b_econconcern b_hardship"		// Post-Treatment Confounders 
local x  = "`controls'"					// Pre-Treatment Confounders 

// Mediators (at most 2 for current code)
local m1 "b_conspiracy"
local m2 "b_polarization"


capture gen ytilde=.

// Both Mediators 
capture program drop deboot1
program define deboot1, rclass
  #delimit ;
  version 17.0;
  args y a m1 x z m2;
  regress `y' `m1' `m2' `x' `z';
  replace ytilde = `y' - _b[`m1']*`m1' - _b[`m2']*`m2';
  reg ytilde `a' `x';
  return scalar deffect = _b[`a'];
  #delimit cr
end

// Mediator 1
capture program drop deboot2
program define deboot2, rclass
  #delimit ;
  version 17.0;
  args y a m1 x z m2;
  regress `y' `m1'  `x' `z';
  replace ytilde = `y' - _b[`m1']*`m1';
  reg ytilde `a' `x';
  return scalar deffect = _b[`a'];
  #delimit cr
end

// Mediator 2
capture program drop deboot3
program define deboot3, rclass
  #delimit ;
  version 17.0;
  args y a m1 x z m2;
  regress `y' `m2'  `x' `z';
  replace ytilde = `y' - _b[`m2']*`m2';
  reg ytilde `a' `x';
  return scalar deffect = _b[`a'];
  #delimit cr
end


capture frame drop acdetab 						// drop if frame exists
frame create acdetab str20(l) bacde lb95 ub95  	// create frame

regress `y' `a' `x'	// Baseline Regression, including pretreatment confounders

frame post acdetab  ///
		("ATE of veteran") ///
		(_b[`a'])        ///
		(r(table)[5,1])  ///
		(r(table)[6,1])

forvalues i = 1/3 {
			
bootstrap deffect=r(deffect), reps(200) seed(12345): deboot`i' `y' `a' `m1' "`x'" "`z'" `m2'

local medlab = cond(`i'==1, "both", ///
			   cond(`i'==2, "conspiracy", "polarization"))

frame post acdetab      ///
	   ("`medlab'")        ///
	   (r(table)[1,1])  ///
	   (r(table)[5,1]) ///
	   (r(table)[6,1])
}	
	
frame acdetab {
		gsort bacde
		capture gen seq=_n
		replace seq = _n
		list seq l bacde lb95 ub95
		
		replace l = "ACDE of veteran | M = conspiracy" if l == "conspiracy"
		replace l = "ACDE of veteran | M = support for Trump" if l=="polarization"
		replace l = "ACDE of veteran | M = both" if l == "both"
		
		twoway  rcap lb95 ub95 seq, horizontal || ///
		       scatter  seq bacde, msymbol(O) mcolor(black) mlabel(l) mlabpos(12) mlabgap(1) mlabsize(medium) ///
			   xline(0, lcolor(red)) ///
			   xlabel(-.05(.05).35, labsize(medium)) scheme(plotplain) yscale(range(.5 4.5)) xtitle(Estimated Effect of Veterancy) ///
			   ylabel(,nolabels notick) ytitle("Mediator Specifications") legend(off) ///
			   xsize(5) ysize(3) 
		graph export fig_15.png, width(3000) replace
}

frame  drop acdetab
	
restore 

********************************************************************************
*
* Figure A.16. Average Controlled Direct Effect of Veteran on Insurrectionist Sentiments by Mediator With Less Restrictive Specifications of IVs and Mediators
*
*********************************************************************************	

// Mediation Analysis, With Continuous Variables 

local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white"


	// The mediators all need to have a 0 that can serve as the baseline 
	foreach var of varlist conspiracy polarization {
		
		gen z`var' = `var' - 1
	}
	  
	// Define Variables Used in the Mediation Analysis, by Category 
	local y  "insurrection_index" 			// Dependent Variable 
	local a  "veteran"						// Treatment 
	local z  "b_econconcern b_hardship"		// Post-Treatment Confounders 
	local x  = "`controls'"					// Pre-Treatment Confounders 

	// Mediators (at most 2 for current code)
	local m1 "zconspiracy"
	local m2 "zpolarization"


	capture gen ytilde=.

	// Both Mediators 
	capture program drop deboot1
	program define deboot1, rclass
	  #delimit ;
	  version 17.0;
	  args y a m1 x z m2;
	  regress `y' `m1' `m2' `x' `z';
	  replace ytilde = `y' - _b[`m1']*`m1' - _b[`m2']*`m2';
	  reg ytilde `a' `x';
	  return scalar deffect = _b[`a'];
	  #delimit cr
	end

	// Mediator 1
	capture program drop deboot2
	program define deboot2, rclass
	  #delimit ;
	  version 17.0;
	  args y a m1 x z m2;
	  regress `y' `m1'  `x' `z';
	  replace ytilde = `y' - _b[`m1']*`m1';
	  reg ytilde `a' `x';
	  return scalar deffect = _b[`a'];
	  #delimit cr
	end

	// Mediator 2
	capture program drop deboot3
	program define deboot3, rclass
	  #delimit ;
	  version 17.0;
	  args y a m1 x z m2;
	  regress `y' `m2'  `x' `z';
	  replace ytilde = `y' - _b[`m2']*`m2';
	  reg ytilde `a' `x';
	  return scalar deffect = _b[`a'];
	  #delimit cr
	end


	capture frame drop acdetab 						// drop if frame exists
	frame create acdetab str20(l) bacde lb95 ub95  	// create frame

	regress `y' `a' `x'	// Baseline Regression, including pretreatment confounders

	frame post acdetab  ///
			("ATE of veteran") ///
			(_b[`a'])        ///
			(r(table)[5,1])  ///
			(r(table)[6,1])

	forvalues i = 1/3 {
				
	bootstrap deffect=r(deffect), reps(200) seed(12345): deboot`i' `y' `a' `m1' "`x'" "`z'" `m2'

	local medlab = cond(`i'==1, "both", ///
				   cond(`i'==2, "conspiracy", "polarization"))

	frame post acdetab      ///
		   ("`medlab'")        ///
		   (r(table)[1,1])  ///
		   (r(table)[5,1]) ///
		   (r(table)[6,1])
	}	
		
	frame acdetab {
			gsort bacde
			capture gen seq=_n
			replace seq = _n
			list seq l bacde lb95 ub95
			
			replace l = "ACDE of veteran | M = conspiracy" if l == "conspiracy"
			replace l = "ACDE of veteran | M = support for Trump" if l=="polarization"
			replace l = "ACDE of veteran | M = both" if l == "both"
			
			twoway  rcap lb95 ub95 seq, horizontal || ///
				   scatter  seq bacde, msymbol(O) mcolor(black) mlabel(l) mlabpos(12) mlabgap(1) mlabsize(medium) ///
				   xline(0, lcolor(red)) ///
				   xlabel(-.05(.05).35, labsize(medium)) scheme(plotplain) yscale(range(.5 4.5)) xtitle(Estimated Effect of Veterancy) ///
				   ylabel(,nolabels notick) ytitle("Mediator Specifications") legend(off) ///
				   xsize(5) ysize(3) 
			graph export fig_16.png, width(3000) replace
	}

	frame  drop acdetab	
	


********************************************************************************
*
* Figure A.17.
* Comparing high Insurrectionist Sentiments (binary DV) across veterans with and without a given pre- or post-military experience, using nonveterans as a benchmark
*
********************************************************************************
estimates clear 


// Change missing to 0 for nonvets (veteran==0) for service and postservice
// treating nonveterans as veterans who do not have the given factor
// for the comparative analyis
   
local servicevars "b_combattrauma b_discharged b_servedgwot b_singletour b_combatdeployment"
local postservicevars "b_unappreciated b_notproud b_reintegration b_disillusioned b_ptsd"


foreach var of varlist `servicevars' `postservicevars' {
	
	replace `var' = 0 if veteran == 0 	// Set all non-veterans to 0	
}

local controls "1.b_male i.age4 1.b_christian 1.b_religious 1.b_childlaw 1.b_republican 1.b_white" 

local counter = 1

foreach iv of varlist `servicevars' `postservicevars' {
	
	logit b_insurrect							/// DV
			i.veteran#i.`iv' 						/// IV
			`controls' 								///
		, robust	
			
	margins i.veteran#i.`iv', atmeans post	// Margins	
	estimates store v1_`counter'			// Store the Estimates for Plot
	
local counter = `counter' + 1
}
		

 * Service Factors (upper panel) ***********************************************

coefplot ///
	(v1_1, keep(1.veteran#1.b*) label("Vet With")  msymbol(d) mcolor(green) ciopts(lcolor(green))) 								///
	(v1_1, keep(1.veteran#0.b*) label("Vet Without") msymbol(s) mcolor(green) ciopts(lcolor(green)))  							///
	(v1_1, keep(0.veteran#0.b*) label("Non-Vet") msymbol(s) mcolor(blue) ciopts(lcolor(blue))), bylabel("Combat" "Trauma") || 	///
	(v1_2, keep(1.veteran#1.b*)) 												///
	(v1_2, keep(1.veteran#0.b*))  												///
	(v1_2, keep(0.veteran#0.b*)), bylabel("Involuntary" "Discharge") || 		///
	(v1_3, keep(1.veteran#1.b*)) 												///
	(v1_3, keep(1.veteran#0.b*))  												///
	(v1_3, keep(0.veteran#0.b*)), bylabel("Service During" "GWoT") 	|| 			///
	(v1_4, keep(1.veteran#1.b*)) 												///
	(v1_4, keep(1.veteran#0.b*))  												///
	(v1_4, keep(0.veteran#0.b*)), bylabel("4 Years or Less" "of Service") || 	///
	(v1_5, keep(1.veteran#1.b*)) 												///
	(v1_5, keep(1.veteran#0.b*))  												///
	(v1_5, keep(0.veteran#0.b*)), bylabel("Combat" "Deployed") 					///
	vertical ///
    msize(huge) mfcolor(white) ///
    mlabel(string(round(@b,0.001), "%9.2f")) /// turn the value into a string and trim the leading zero
    mlabposition(0) ///
    grid(within) ///
	byopts(compact row(1) legend(off)) ///
    drop(_cons) ///
	rename(												///
		1.veteran#1.b_combattrauma     	= "Vet With" 	///
		1.veteran#0.b_combattrauma     	= "Vet Without" ///
		0.veteran#0.b_combattrauma     	= "Non-Vet" 	///
		1.veteran#1.b_discharged 		= "Vet With" 	///
		1.veteran#0.b_discharged 		= "Vet Without" ///
		0.veteran#0.b_discharged 		= "Non-Vet" 	///
		1.veteran#1.b_servedgwot 		= "Vet With" 	///
		1.veteran#0.b_servedgwot 		= "Vet Without" ///
		0.veteran#0.b_servedgwot 		= "Non-Vet" 	///
		1.veteran#1.b_singletour 		= "Vet With" 	///
		1.veteran#0.b_singletour 		= "Vet Without" ///
		0.veteran#0.b_singletour 		= "Non-Vet" 	///
		1.veteran#1.b_combatdeployment 	= "Vet With" 	///
		1.veteran#0.b_combatdeployment 	= "Vet Without"	///
		0.veteran#0.b_combatdeployment	= "Non-Vet" 	///
		) 												///
	title("")  											///
	ytitle("Predicted Proportion") 					///
	xlabel("") 											///
	yline(.0527705, lcolor(blue)) 						/// put line at non-veteran mean value 
	ylab(0(0.05).2) 									///
	legend(off) 										///
	xsize(6.5) ysize(3) scheme(white_tableau) 			///
	name(int_1, replace)
	

* Post-Service Factors (lower panel) *******************************************
	
coefplot ///
	(v1_6, keep(1.veteran#1.b*) label("Vet With")  msymbol(d) mcolor(green) ciopts(lcolor(green))) 										///
	(v1_6, keep(1.veteran#0.b*) label("Vet Without") msymbol(s) mcolor(green) ciopts(lcolor(green)))  									///
	(v1_6, keep(0.veteran#0.b*) label("Non-Vet") msymbol(s) mcolor(blue) ciopts(lcolor(blue))), bylabel("Service" "Unappreciated") || 	///
	(v1_7, keep(1.veteran#1.b*)) 												///
	(v1_7, keep(1.veteran#0.b*))  												///
	(v1_7, keep(0.veteran#0.b*)), bylabel("No Pride in" "Service") || 			///
	(v1_8, keep(1.veteran#1.b*)) 												///
	(v1_8, keep(1.veteran#0.b*))  												///
	(v1_8, keep(0.veteran#0.b*)), bylabel("Reintegration" "Trouble") || 		///
	(v1_9, keep(1.veteran#1.b*)) 												///
	(v1_9, keep(1.veteran#0.b*))  												///
	(v1_9, keep(0.veteran#0.b*)), bylabel("Disillusioned" "with Service") || 	///
	(v1_10, keep(1.veteran#1.b*)) 												///
	(v1_10, keep(1.veteran#0.b*)) 	 											///
	(v1_10, keep(0.veteran#0.b*)), bylabel("Service-Related" "PTSD") 			///
	vertical ///
    msize(huge) mfcolor(white) ///
    mlabel(string(round(@b,0.001), "%9.2f")) /// turn the value into a string and trim the leading zero
    mlabposition(0) ///
    grid(within) ///
	byopts(compact row(1) legend(off)) ///
    drop(_cons) ///
	rename(	///
		1.veteran#1.b_notproud     	= "Vet With" ///
		1.veteran#0.b_notproud     	= "Vet Without" ///
		0.veteran#0.b_notproud     	= "Non-Vet" ///
		1.veteran#1.b_reintegration = "Vet With" ///
		1.veteran#0.b_reintegration = "Vet Without" ///
		0.veteran#0.b_reintegration = "Non-Vet" ///
		1.veteran#1.b_disillusioned = "Vet With" ///
		1.veteran#0.b_disillusioned = "Vet Without" ///
		0.veteran#0.b_disillusioned = "Non-Vet" ///
		1.veteran#1.b_unappreciated = "Vet With" ///
		1.veteran#0.b_unappreciated = "Vet Without" ///
		0.veteran#0.b_unappreciated = "Non-Vet" ///
		1.veteran#1.b_ptsd          = "Vet With" ///
		1.veteran#0.b_ptsd          = "Vet Without" ///
		0.veteran#0.b_ptsd          = "Non-Vet") ///
	title("")  ///
	ytitle("Predicted Proportion") ///
	xlabel("") ///
	yline(.0527705, lcolor(blue)) /// put line at non-veteran mean value 
	xsize(6.5) ysize(3) scheme(white_tableau) ylab(0(.05).2) ///
	name(int_2, replace) 

* Combine top and bottom *******************************************************	

graph combine int_1 int_2, rows(2) ycommon xcommon
graph export fig_a17.png, replace width(3000)


********************************************************************************
*
* Figure A.7
* Table A7. ACDE of Veteran Across Alternative Specifications of Mediation Analysis
* 
******************************************************************************** 	 

putexcel set tab_a7.xlsx, replace 

putexcel A2 = "#"
putexcel B2 = "Pretreatment Confounders (X)"
putexcel C2 = "Intermediate Variables (Z)"

putexcel D2 = "ATE of Veteran"

putexcel E1:F1 = "Both", merge hcenter
putexcel E2 = "ACDE"
putexcel F2 = "Sig", italic

putexcel G1:H1 = "Conspiracy", merge hcenter
putexcel G2 = "ACDE"
putexcel H2 = "Sig", italic

putexcel I1:J1 = "Support for Trump", merge hcenter
putexcel I2 = "ACDE"
putexcel J2 = "Sig", italic

* Define Variables Used in the Mediation Analysis, by Category *****************

local y  "insurrection_index" 			// Dependent Variable 
local a  "veteran"						// Treatment 
local z  "b_econconcern b_hardship"		// Post-Treatment Confounders 

// Pre-Treatment Confounders 
local x  "b_republican b_childlaw b_white b_under35 b_male b_christian b_religious"

// Mediators (at most 3 for current code)
local m1 "b_conspiracy"
local m2 "b_polarization"


* Mediation Analysis Implementation ********************************************

capture gen ytilde=.

// Both Mediators 
capture program drop deboot1
program define deboot1, rclass
  #delimit ;
  version 17.0;
  args y a m1 x z m2;
  regress `y' `m1' `m2' `x' `z';
  replace ytilde = `y' - _b[`m1']*`m1' - _b[`m2']*`m2';
  reg ytilde `a' `x';
  return scalar deffect = _b[`a'];
  #delimit cr
end

// Mediator 1
capture program drop deboot2
program define deboot2, rclass
  #delimit ;
  version 17.0;
  args y a m1 x z m2;
  regress `y' `m1'  `x' `z';
  replace ytilde = `y' - _b[`m1']*`m1';
  reg ytilde `a' `x';
  return scalar deffect = _b[`a'];
  #delimit cr
end

// Mediator 2
capture program drop deboot3
program define deboot3, rclass
  #delimit ;
  version 17.0;
  args y a m1 x z m2;
  regress `y' `m2'  `x' `z';
  replace ytilde = `y' - _b[`m2']*`m2';
  reg ytilde `a' `x';
  return scalar deffect = _b[`a'];
  #delimit cr
end



local currentrow = 3

*drop b_educ1 // Otherwise the bootstrap has trouble with the estimations

foreach variant of numlist 1/16 {

	if `variant' == 1  	local z "b_econconcern b_hardship"
	if `variant' == 1	local x "b_republican b_childlaw b_white b_under35 b_male b_christian b_religious"
	if `variant' == 2 	local z "b_econconcern"
	if `variant' == 2	local x "b_republican b_childlaw b_white b_under35 b_male b_christian b_religious"
	if `variant' == 3  	local z "b_hardship"
	if `variant' == 3	local x "b_republican b_childlaw b_white b_under35 b_male b_christian b_religious"
	if `variant' == 4  	local z "b_econconcern b_hardship"
	if `variant' == 4	local x "b_childlaw b_white b_under35 b_male b_christian b_religious"
	if `variant' == 5  	local z "b_econconcern b_hardship"
	if `variant' == 5	local x "b_republican b_white b_under35 b_male b_christian b_religious"
	if `variant' == 6  	local z "b_econconcern b_hardship"
	if `variant' == 6	local x "b_republican b_childlaw b_under35 b_male b_christian b_religious"
	if `variant' == 7  	local z "b_econconcern b_hardship"
	if `variant' == 7	local x "b_republican b_childlaw b_white b_male b_christian b_religious"
	if `variant' == 8  	local z "b_econconcern b_hardship"
	if `variant' == 8	local x "b_republican b_childlaw b_white b_under35 b_christian b_religious"
	if `variant' == 9  	local z "b_econconcern b_hardship"
	if `variant' == 9	local x "b_republican b_childlaw b_white b_under35 b_male b_religious"
	if `variant' == 10  local z "b_econconcern b_hardship"
	if `variant' == 10	local x "b_republican b_childlaw b_white b_under35 b_male b_christian"
	if `variant' == 11  local z "b_econconcern b_hardship b_republican" 
	if `variant' == 11	local x "b_childlaw b_white b_under35 b_male b_christian b_religious"
	if `variant' == 12  local z "b_econconcern b_hardship b_christian" 
	if `variant' == 12	local x "b_republican b_childlaw b_white b_under35 b_male b_religious"
	if `variant' == 13  local z "b_econconcern b_hardship b_religious" 
	if `variant' == 13	local x "b_republican b_childlaw b_white b_under35 b_male b_christian"
	if `variant' == 14  local z "b_econconcern b_hardship b_republican b_christian" 
	if `variant' == 14	local x "b_childlaw b_white b_under35 b_male  b_religious"
	if `variant' == 15  local z "b_econconcern b_hardship b_christian b_religious" 
	if `variant' == 15	local x "b_republican b_childlaw b_white b_under35 b_male"
	if `variant' == 16  local z "b_econconcern b_hardship b_republican b_christian b_religious" 
	if `variant' == 16	local x "b_childlaw b_white b_under35 b_male"
	
	if `variant' == 1   local xlab = "G C W Y M X R"
	if `variant' == 2 local xlab = "G C W Y M X R"
	if `variant' == 3 local xlab = "G C W Y M X R"
	if `variant' == 4 local xlab = "C W Y M X R"
	if `variant' == 5 local xlab = "G W Y M X R"
	if `variant' == 6 local xlab = "G C Y M X R"
	if `variant' == 7 local xlab = "G C W M X R"
	if `variant' == 8 local xlab = "G C W Y X R"
	if `variant' == 9 local xlab = "G C W Y M R"
	if `variant' == 10 local xlab = "G C W Y M X"
	if `variant' == 11 local xlab = "C W Y M X R"
	if `variant' == 12 local xlab = "G C W Y M R"
	if `variant' == 13 local xlab = "G C W Y M X"
	if `variant' == 14 local xlab = "C W Y M R"
	if `variant' == 15 local xlab = "G C W Y M"
	if `variant' == 16 local xlab = "C W Y M"
	if `variant' == 1   local zlab "E H"
	if `variant' == 2 local zlab "E"
	if `variant' == 3 local zlab "H"
	if `variant' == 4 local zlab "E H"
	if `variant' == 5 local zlab "E H"
	if `variant' == 6 local zlab "E H"
	if `variant' == 7 local zlab "E H"
	if `variant' == 8 local zlab "E H"
	if `variant' == 9 local zlab "E H"
	if `variant' == 10 local zlab "E H"
	if `variant' == 11 local zlab "E H G"
	if `variant' == 12 local zlab "E H X"
	if `variant' == 13 local zlab "E H R"
	if `variant' == 14 local zlab "E H G X"
	if `variant' == 15 local zlab "E H X R"
	if `variant' == 16 local zlab "E H G X R"


	
	regress `y' `a' `x'	// Baseline Regression, including pretreatment confounders

	putexcel A`currentrow' = `variant'
	putexcel B`currentrow' = "`xlab'"
	putexcel C`currentrow' = "`zlab'"
	putexcel D`currentrow' = _b[`a']
	

		forvalues i = 1/3 {
					
			bootstrap deffect=r(deffect), reps(200) seed(12345): deboot`i' `y' `a' `m1' "`x'" "`z'" `m2'
			
		//	local medlab = cond(`i'==1, "both", ///
		//	   cond(`i'==2, "conspiracy", "polarization"))

						local col1 = "E"
						local col2 = "F"
			if `i' == 2 local col1 = "G"
			if `i' == 2 local col2 = "H"
			if `i' == 3 local col1 = "I"
			if `i' == 3 local col2 = "J"
			
			local acde = r(table)[1,1]
			local sig = r(table)[4,1]

				if `sig' < 0.05 	local star = "*"
				if `sig' < 0.01 	local star = "**"
				if `sig' < 0.001 	local star = "***"
			
			putexcel `col1'`currentrow' = `acde'
			putexcel `col2'`currentrow' = "`star'"		
		}
	
	local currentrow = `currentrow' + 1		
	
}

	local nextrow = `currentrow' + 1

putexcel A`currentrow':J`nextrow' = "C - Childhood Trouble with the Law, E - Economic Concern, H - Recent Economic Hardship, G - Republican, M - Male, R - Religious, W - White, X - Christian, Y - Under 35", merge left txtwrap

	local nextrow = `nextrow' + 1

putexcel A`nextrow':J`nextrow' = "* p < 0.05, ** p < .01, *** p < .001", merge left txtwrap

putexcel save


// End of dofile
